﻿using System;
using AGLibrary;

namespace ErfExpIntCalculator
{
    public class fsErfExpIntCalculator
    {
        /*
         * This class was automatically generated from the Maple module using CodeGeneration[CSharp] procedure
         * and then slightly tuned.
         * What's happening here is a terrible mystery! :-)
         * See the attached Maple doc (if I'll have a time and a possibility to attach it!).
         * Also see the comments in fsSpecialFunctions.cs
         * ---Tyshkevich Dmitry---.
         */
        public static double ErfExpInt(double a, double b, double x)
        {  
            if (x <= -xT) 
                return 0.0e0;
            if (xT < x)
                return ErfExpInt(a, b, xT);
            if (a == 0.0e0)
                return normaldistr.erf(b) * (1 + normaldistr.erf(x));
            if (a < 0.0e0)
                return -ErfExpInt(-a, -b, x);
            if (0.1e1 < a)
                return normaldistr.erf(a * x + b) * normaldistr.erf(x) - 1 - ErfExpInt(0.1e1 / a, -b / a, a * x + b);
            return ErfExpIntSpecial(a, b, x);
            
        }

        public static double ErfcExpInt(double a, double b, double x)
        {
            return 1 + normaldistr.erf(x) - ErfExpInt(a, b, x);
        }

        private static readonly double PiSqrt = 0.177245385090551602729816748334e1;

        private static readonly double epsHalf = 0.5000000000e-9;

        private static readonly double epsQuarter = 0.2500000000e-9;

        private static readonly double xT = 0.4397571017214433973343322e1;

        private static readonly int nIter = 15;

        private static readonly int nBounds = 40;

        private static readonly double[] bounds = new[]{
  0.5e0, 0.7375e0, 0.975e0, 0.12125e1, 0.145e1, 0.16875e1, 0.1925e1, 0.21625e1, 0.24e1, 0.26375e1, 0.2875e1, 
  0.31125e1, 0.335e1, 0.35875e1, 0.3825e1, 0.40625e1, 0.43e1, 0.45375e1, 0.4775e1, 0.50125e1, 0.525e1, 0.54875e1, 
  0.5725e1, 0.59625e1, 0.62e1, 0.64375e1, 0.6675e1, 0.69125e1, 0.715e1, 0.73875e1, 0.7625e1, 0.78625e1, 0.81e1, 
  0.83375e1, 0.8575e1, 0.88125e1, 0.905e1, 0.92875e1, 0.9525e1, 0.97625e1, 0.100e2};

        private static readonly double[] shifts = new[]{
  0.618750000014551915228367e0, 0.856249999997089616954325e0, 0.109374999997962731868027e1, 
  0.133124999999781721271575e1, 0.156875000001600710675121e1, 0.180624999999854480847715e1, 
  0.204374999998108251020312e1, 0.228125000003492459654809e1, 0.251875000001746229827404e1, 
  0.275625e1, 0.299375000000036379788071e1, 0.323125000000072759576141e1, 0.346875000001891748979688e1, 
  0.370625000000145519152285e1, 0.394374999998399289324879e1, 0.418125000000218278728425e1, 
  0.441875000002037268131972e1, 0.465625000000291038304568e1, 0.489374999998544808477163e1, 
  0.513125000000363797880709e1, 0.536874999995052348822354e1, 0.560625000000436557456851e1, 
  0.584374999998690327629447e1, 0.608125000004074536263943e1, 0.631874999998763087205588e1, 
  0.655625000000582076609136e1, 0.679374999998835846781731e1, 0.703124999997089616954326e1, 
  0.726874999995343387126923e1, 0.750625000000727595761418e1, 0.774374999991850927472115e1, 
  0.798124999997235136106611e1, 0.821875000002619344741106e1, 0.845625000000873114913701e1, 
  0.869374999999126885086299e1, 0.893125000011641532182694e1, 0.916875000009895302355289e1, 
  0.940624999993888195604085e1, 0.964375000006402842700482e1, 0.979999999992869561538101e1};

        private static readonly double[] coeffsErf = new[]{
  -0.333333333333333333333333e0, 0.1e0, -0.238095238095238095238095e-1, 0.462962962962962962962963e-2, 
  -0.757575757575757575757576e-3, 0.106837606837606837606838e-3, -0.132275132275132275132275e-4, 
  0.145891690009337068160598e-5, -0.145038522231504687645039e-6, 0.131225329638028050726463e-7, 
  -0.108922210371485733804574e-8};

        private static readonly double[] coeffsErfc = new[]{
 0.279763689203807902590448e0, -0.217982018149901814442494e0, 0.144887315470384083801509e0, -0.855553278003287423144249e-1, 0.459749781963432104545924e-1, -0.228433240162709044609932e-1, 0.106135571536455474504835e-1, -0.465033872206882830464019e-2, 0.193404001731822800336980e-2, -0.767478102576552961543818e-3, 0.291832588362789377287746e-3, -0.106710306971151306198947e-3, 0.376342640642080067866884e-4, -0.128344775551091464485034e-4, 0.424184727284258099165650e-5, 0.235150109394083386202506e0, -0.161495021211757238879551e0, 0.968699974819862606650986e-1, -0.523667239120590745645659e-1, 0.260154950662190722750939e-1, -0.120363825046725479344295e-1, 0.523644751554204805915213e-2, -0.215790694843627437380231e-2, 0.847184922743831938900050e-3, -0.318334412951553823377704e-3, 0.114922216307249727382316e-3, -0.399876845905263025664474e-4, 0.134471269317238940638566e-4, -0.438055113363956245314393e-5, 0.138518285645882572011403e-5, 0.201632118166514805327476e0, -0.123119325066720727819762e0, 0.669703563772972898672616e-1, -0.332469985202773298161581e-1, 0.153032258732105602932743e-1, -0.660363808870709874409653e-2, 0.269349890460771744994932e-2, 
 -0.104503533195247386246820e-2, 0.387622877574061798639728e-3, -0.138016179920842979394919e-3, 0.473335361646699804631983e-4, -0.156809317671394319059619e-4, 0.503041949790163397729984e-5, -0.156598621842161427392112e-5, 0.473946014629587440373476e-6, 0.175767574424705900794699e0, -0.962084166427441463333244e-1, 0.476901197692627498604542e-1, -0.218139631333449728249777e-1, 0.932514067402240242288066e-2, -0.375994784442808661954893e-2, 0.143990336871232264664505e-2, -0.526593281383983482145687e-3, 0.184719015717114948320699e-3, -0.623746870421650922537191e-4, 0.203365427178848345497178e-4, -0.641848446734029580193805e-5, 0.196532254728973261986751e-5, -0.584945972109361593510452e-6, 0.169516172756145843176127e-6, 0.155342931858314743410168e0, -0.768011348373206098487691e-1, 0.348611515810386564452705e-1, -0.147418021960052330387184e-1, 0.586747469291014854664445e-2, -0.221488050856436457328993e-2, 0.797626965019582932864956e-3, -0.275315202048980446357480e-3, 0.914315604533349580528116e-4, -0.293070981320953405677930e-4, 0.909121004925856790121907e-5, -0.273550224694010709185540e-5, 0.799981817231186159595428e-6, -0.227773965752852161337859e-6, 
 0.632373442685973990075179e-7, 0.138888687790764631842534e0, -0.624541989040232537543967e-1, 0.260807910204636184979627e-1, -0.102305134155660077850431e-1, 0.380096308180512070716896e-2, -0.134600953962478531207854e-2, 0.456577783622729224414301e-3, -0.148947405128564812312443e-3, 0.468853832745387955706811e-4, -0.142801514663335149138110e-4, 0.421837193935990310102612e-5, -0.121103948120113407710649e-5, 0.338488645321992395264238e-6, -0.922529793114557575258214e-7, 0.245509575504006254537284e-7, 0.125401034576743618719275e0, -0.516128547200612949550253e-1, 0.199172627435944685002289e-1, -0.727129932547847696019234e-2, 0.252827237364370055733109e-2, -0.841657064753459808326050e-3, 0.269378582523739951194803e-3, -0.831755962119602111258864e-4, 0.248471144392750156062058e-4, -0.719873468154425674244479e-5, 0.202694008816737504237646e-5, -0.555668341880557725050606e-6, 0.148548818877887854086856e-6, -0.387802602925680859362187e-7, 0.989880880889018260994593e-8, 0.114176482988066348064658e0, -0.432593799067285590257363e-1, 0.154910225743310770306609e-1, -0.528032310566134474814446e-2, 0.172264274468127761164783e-2, -0.540217737721021548480547e-3, 
 0.163423676820817283499025e-3, -0.478307071370601835367301e-4, 0.135774690447802467795754e-4, -0.374602352536472652834170e-5, 0.100637057490607914557377e-5, -0.263680118998243182080001e-6, 0.674750507784290471347143e-7, -0.168850244816296750792458e-7, 0.413658405042844300439286e-8, 0.104710533065329465531565e0, -0.367102732275018248305739e-1, 0.122465323729198594802995e-1, -0.390954654199631967863280e-2, 0.119968101009232138365815e-2, -0.355139999135538749051760e-3, 0.101724045750322172201738e-3, -0.282635882478562999051125e-4, 0.763378322101222498410164e-5, -0.200799927934397226302790e-5, 0.515227005955555651382232e-6, -0.129140228719681229813069e-6, 0.316591759797629508793424e-7, -0.759979684803901494712947e-8, 0.178816965906717441962472e-8, 0.966338157090799532204839e-1, -0.314956744514523544236022e-1, 0.982386300226269522964336e-2, -0.294576803432326762095069e-2, 0.852294928809538613927733e-3, -0.238652054718798025486066e-3, 0.648367343547189884254980e-4, -0.171273730193346149095463e-4, 0.440735312427421279847659e-5, -0.110657910651242462094352e-5, 0.271468891626292094422575e-6, -0.651532676973777073633971e-7, 0.153150329502270617367515e-7, 
 -0.352941667297326651436110e-8, 0.798154043007415390520042e-9, 0.896707996925723396816592e-1, -0.272856703884173967294115e-1, 0.798432396722158254719476e-2, -0.225506700771997224137435e-2, 0.616608556440583641435521e-3, -0.163638056684256007801919e-3, 0.422390414872996572442479e-4, -0.106242646051312410745962e-4, 0.260816233272243904844875e-5, -0.625795252061499683425447e-6, 0.146937559368885115693927e-6, -0.338001695137907923991406e-7, 0.762471693491375842054800e-8, -0.168825742457463699093242e-8, 0.367213759774200447824737e-9, 0.836124541212926360578311e-1, -0.238440987887938524387347e-1, 0.656620990996883767828850e-2, -0.175135534475389497402967e-2, 0.453571476282495242779158e-3, -0.114301004597487053207733e-3, 0.280787851851606334307356e-4, -0.673469424711534795199380e-5, 0.157932609550691560284123e-5, -0.362554846882876039169761e-6, 0.815641497584507208330573e-7, -0.180001249874216486519289e-7, 0.390020759709145311661892e-8, -0.830396793636643319597591e-9, 0.173855425017520361387377e-9, 0.782977554405582025274504e-1, -0.209989051759210887901292e-1, 0.545780311129315915069026e-2, -0.137810042198578638515504e-2, 0.338758636786096306178606e-3, 
 -0.812125598314915052195364e-4, 0.190175234559277801732602e-4, -0.435586436622154807945650e-5, 0.977029726159474185003204e-6, -0.214842779539835270513730e-6, 0.463587674707024850076142e-7, -0.982469181489432828288530e-8, 0.204656125972507319833835e-8, -0.419335752849811623327024e-9, 0.845700522541319667177350e-10, 0.736009276529220128989519e-1, -0.186227073201087536307160e-1, 0.458051864835246187545377e-2, -0.109744005194625023036244e-2, 0.256565729092018363499654e-3, -0.586173267282217756330814e-4, 0.131050874843962016727347e-4, -0.287045606651402898777216e-5, 0.616614913488914402626699e-6, -0.130028231666001906680158e-6, 0.269395563877477237583590e-7, -0.548790930904967424281974e-8, 0.109999874363933715663981e-8, -0.217082912094800687005601e-9, 0.422050277085792011206807e-10, 0.694224769539922269856357e-1, -0.166197965746566855985837e-1, 0.387815421488441465418955e-2, -0.883550590480090048905034e-3, 0.196825788976369378131666e-3, -0.429275531557478768284941e-4, 0.917675057903690123272922e-5, -0.192478374468150964450747e-5, 0.396471164321591812499999e-6, -0.802667969225110984480320e-7, 0.159837970043832690628208e-7, -0.313285416365530040532703e-8, 
 0.604767237120756009712524e-9, -0.115046670092847122892678e-9, 0.215788465980360116124780e-10, 0.656827862961665866107025e-1, -0.149172831435656847430908e-1, 0.330989615680446879413408e-2, -0.718519886656621887961226e-3, 0.152792443729530166255026e-3, -0.318625913722131057925377e-4, 0.652232799187187404314068e-5, -0.131174501223506557235708e-5, 0.259398535089484843082689e-6, -0.504744181884806965349463e-7, 0.967047484753139966060660e-8, -0.182540827905110943048955e-8, 0.339664422127479609841129e-9, -0.623363709637226692909325e-10, 0.112886382580892026769294e-10, 0.623174679136732228796694e-1, -0.134589608525899811509582e-1, 0.284568465579165585903681e-2, -0.589727843992569403244977e-3, 0.119912377413350632234592e-3, -0.239460091976175628330945e-4, 0.470031651673624242565690e-5, -0.907567342513016569086443e-6, 0.172500828347621333568879e-6, -0.322954004291890260494408e-7, 0.595910546880271304799286e-8, -0.108429151005996714136611e-8, 0.194648734618072496009384e-9, -0.344903719679723799019591e-10, 0.603491461346992525454058e-11, 0.592739430645250250808234e-1, -0.122009887470052665672892e-1, 0.246308922937905759290774e-2, -0.488153000590205745060366e-3, 
 0.950634166659817717377381e-4, -0.182055852808165996499976e-4, 0.343122010629863689229177e-5, -0.636847650527969066176431e-6, 0.116474558521956489981197e-6, -0.210028851502150011288731e-7, 0.373597478558146108682318e-8, -0.655864111058153630948035e-9, 0.113684592218957837556762e-9, -0.194646513931584122850041e-10, 0.329318701806410466162400e-11, 0.565088916228594079134584e-1, -0.111088067672038092386732e-1, 0.214516853689576139844340e-2, -0.407258805032696619226338e-3, 0.760728880602422568969490e-4, -0.139908422041660588786794e-4, 0.253506799295790211834979e-5, -0.452815235072828757200653e-6, 0.797758609077226330992465e-7, -0.138693584239183044219786e-7, 0.238053749333423449075807e-8, -0.403564191089456579773348e-9, 0.675992088340939366524334e-10, -0.111923950248501583502571e-10, 0.183234658558397976227330e-11, 0.539863269947951357222979e-1, -0.101549027210514117821879e-1, 0.187898245651821672806016e-2, -0.342249297290310973660998e-3, 0.614078847448091183063961e-4, -0.108600341025479467653516e-4, 0.189411152207152863367069e-5, -0.325964128022375255092348e-6, 0.553770253966768813495798e-7, -0.929128109556134669076124e-8, 0.154022774093710224058878e-8, 
 -0.252361354417683732333008e-9, 0.408830915407822297129674e-10, -0.655076835087803841456027e-11, 0.103849456184662784972537e-11, 0.516761265154231795375548e-1, -0.931717502230070501015018e-2, 0.165454318899809784993182e-2, -0.289564144122580891452810e-3, 0.499728565411098376513639e-4, -0.850894697267537669053026e-5, 0.143014905794263035267388e-5, -0.237381225793011338769578e-6, 0.389271538093809856600323e-7, -0.630912629619377981567228e-8, 0.101100631467260344696338e-8, -0.160233833275712969821471e-9, 0.251251538365705136655481e-10, -0.389894854517093685884486e-11, 0.598953462577493524091114e-12, 0.495528976964549205009114e-1, -0.857771801198721449394718e-2, 0.146406619845681232562163e-2, -0.246531206223128674250264e-3, 0.409753245924374538747230e-4, -0.672531600367637789748513e-5, 0.109050711713799863008480e-5, -0.174760150934699265269770e-6, 0.276895085719587288047940e-7, -0.433907640607793351964128e-8, 0.672712235970297473052255e-9, -0.103215159210994934257112e-9, 0.156770410792668046837033e-10, -0.235780746788087960216582e-11, 0.351226176612964748409551e-12, 0.475950949776343225563482e-1, -0.792191082391204646964188e-2, 0.130142874855391912709378e-2, 
 -0.211124318601645067977477e-3, 0.338355198462225551177318e-4, -0.535919891561440954185889e-5, 0.839233577299422272735550e-6, -0.129979350565815195207735e-6, 0.199166899107363125082639e-7, -0.302026549059476415314971e-8, 0.453402656863866830148770e-9, -0.673988486671347158318478e-10, 0.992343838589768525025566e-11, -0.144750089512645337325368e-11, 0.209229312846763415631692e-12, 0.457843252971774362586463e-1, -0.733772686112270237509080e-2, 0.116177402148687944901196e-2, -0.181792317088199418959275e-3, 0.281247614955883203756551e-4, -0.430344383398176874446653e-5, 0.651481067214555524192516e-6, -0.976070296533949716475544e-7, 0.144770821353060190497573e-7, -0.212628332684575178626627e-8, 0.309324316928314985538820e-9, -0.445826891035086097695317e-10, 0.636763901880262619010601e-11, -0.901459104670334719933628e-12, 0.126520130897063046745187e-12, 0.441047974505794118306810e-1, -0.681520540269430977076922e-2, 0.104121857203058839588453e-2, -0.157336941120624493157588e-3, 0.235229032168781560401890e-4, -0.348063813358066905822038e-5, 0.509873509834097313155377e-6, -0.739642535902303688073711e-7, 0.106279726951985464968587e-7, -0.151305605477724533853489e-8, 
 0.213469949534180027391776e-9, -0.298532342259426992849635e-10, 0.413913695354515017737633e-11, -0.569086508850304967140277e-12, 0.776030859866351500234933e-13, 0.425428819408561054449225e-1, -0.634604359347839616974949e-2, 0.936633962061226730935816e-3, -0.136824678822602175732373e-3, 0.197885966995751499207118e-4, -0.283427647523086351149538e-5, 0.402123698396130400382871e-6, -0.565294170840713850152956e-7, 0.787567846296072831304930e-8, -0.108766694188792310088907e-8, 0.148932423311756107850131e-9, -0.202234050996159769398137e-10, 0.272378674685796939289343e-11, -0.363935063906858345529170e-12, 0.482482116833264594538102e-13, 0.410867565160723338620743e-1, -0.592327870598243335955849e-2, 0.845482220472241130958900e-3, -0.119522457574092334793927e-3, 0.167382782198147614787699e-4, -0.232271197136450540626629e-5, 0.319451110905734405614480e-6, -0.435545600413446504240501e-7, 0.588808072365721082477175e-8, -0.789425019133218566390345e-9, 0.104984912907564169278484e-9, -0.138515928365736626424464e-10, 0.181344197012665478252216e-11, -0.235618639607612974870485e-12, 0.303868316913896074025890e-13, 0.397261186802412061974146e-1, -0.554103871070625327303776e-2, 
 0.765690751816049230207311e-3, -0.104850269937072001803904e-3, 0.142311615943048717352186e-4, -0.191496619036125281033561e-5, 0.255518392950859824939367e-6, -0.338149879340430426263298e-7, 0.443919040582776746621186e-8, -0.578206934276910144737107e-9, 0.747345946034528551021109e-10, -0.958714806574008245298996e-11, 0.122082637937808415336261e-11, -0.154340361040542460365035e-12, 0.193743851956446621509970e-13, 0.384519509817188009789632e-1, -0.519434499153280240490307e-2, 0.695556434563128381580405e-3, -0.923459521559392126171015e-4, 0.121584130070949970443111e-4, -0.158779541754925175230393e-5, 0.205708241547130738179895e-6, -0.264438861579628338025221e-7, 0.337356130760496459043383e-8, -0.427180679021938041716426e-9, 0.536983641313281382296617e-10, -0.670194476563269118375618e-11, 0.830600370305313868787873e-12, -0.102233560708732057830154e-12, 0.124985938971819150865662e-13, 0.372563281850098346883087e-1, -0.487895519029758611587668e-2, 0.633671511890878321550410e-3, -0.816387684701812700439739e-4, 0.104352680861309569410926e-4, -0.132361548248576432873479e-5, 0.166626360231602217487811e-6, -0.208218084860593056830265e-7, 0.258316509576833532962510e-8, 
 -0.318205774671324409967160e-9, 0.389266135083865565643239e-10, -0.472961549601065502319688e-11, 0.570822788333991071422017e-12, -0.684425869179535065469579e-13, 0.815365779661289786357902e-14, 0.361322579193383743121490e-1, -0.459123717532849745289668e-2, 0.578865891813115720644028e-3, -0.724294354776295316814370e-4, 0.899523977713967941984933e-5, -0.110901963020481678748920e-5, 0.135756259101314006827212e-6, -0.165020180155649426644860e-7, 0.199218909421970867083346e-8, -0.238889795840624682602428e-9, 0.284572593793194263156558e-10, -0.336798136754916298887466e-11, 0.396075583102410108475733e-12, -0.462878442985689170197740e-13, 0.537629665793206734681040e-14, 0.350735482498683699378763e-1, -0.432806732405437172981698e-2, 0.530161902603277282090787e-3, -0.644748977364401869491282e-4, 0.778582604794197446593478e-5, -0.933710186525159068209011e-6, 0.111217144357845125536570e-6, -0.131595197373437392253806e-7, 0.154693153974848038669845e-8, -0.180682807454584276273078e-9, 0.209713864329678091691199e-10, -0.241907884591329583932064e-11, 0.277352235970374002617374e-12, -0.316094249426567063331282e-13, 0.358135774465285201962934e-14, 0.340746970537947237234627e-1, 
 -0.408674793661993844453913e-2, 0.486738575638957029010607e-3, -0.575766380985002330931647e-4, 0.676530321833315804634577e-5, -0.789721737346911453176454e-6, 0.915925158699499785039562e-7, -0.105559188540622696533035e-7, 0.120901402880422389591632e-8, -0.137629990536950532435368e-9, 0.155735167574872024848356e-10, -0.175184611604540559970530e-11, 0.195921936942421102370054e-12, -0.217865645442523363417404e-13, 0.240908620013967689191278e-14, 0.331307991363432696861126e-1, -0.386493979556219134155502e-2, 0.447903267969726010891168e-3, -0.515716384308902355181739e-4, 0.590031181018629473067053e-5, -0.670851559711016515390163e-6, 0.758077364561288415377975e-7, -0.851495992393449261303297e-8, 0.950776177777797424069170e-9, -0.105546432896505265077050e-9, 0.116498370787001216599276e-10, -0.127863666021070699998426e-11, 0.139561000837872588723359e-12, -0.151498361987527216449162e-13, 0.163574206014779668137191e-14, 0.322374678791628749181758e-1, -0.366060679654612976512463e-2, 0.413068979267173818410259e-3, -0.463253427880163509971791e-4, 0.516402569275463857093409e-5, -0.572238691719665432962913e-6, 0.630418353481593109590838e-7, -0.690534580621255964941126e-8, 
 0.752120816802125193982160e-9, -0.814656611776205630228600e-10, 0.877574958950675310422299e-11, -0.940271128160485545288892e-12, 0.100211278491651320631700e-12, -0.106245114115180236042165e-13, 0.112063284512319678076529e-14, 0.313907688487470210340656e-1, -0.347197025470177631655214e-2, 0.381736115344023967424054e-3, -0.417261443107248001167086e-4, 0.453475386270535735761388e-5, -0.490050493183275771629678e-6, 0.526634561960645310395967e-7, -0.562856513781739077963517e-8, 0.598332958452697222427653e-9, -0.632675290585355867441576e-10, 0.665497120001102046527712e-11, -0.696421822477545631352637e-12, 0.725089991953261996587294e-13, -0.751166579883254739774870e-14, 0.774347519808131271494220e-15, 0.305871632994094591103796e-1, -0.329747101085695818229976e-2, 0.353477745704746980616613e-3, -0.376810395334049865006033e-4, 0.399486549499071300973293e-5, -0.421247559590281828784327e-6, 0.441839784115254957146635e-7, -0.461019677507032655051070e-8, 0.478558706554144587272576e-9, -0.494247962728435228561927e-10, 0.507902336999683526659502e-11, -0.519364136332537289943274e-12, 0.528506041765450569771389e-13, -0.535233332917004850841225e-14, 0.539485330544100806924498e-15, 
 0.298234598952170088166498e-1, -0.313573785819943878852925e-2, 0.327927622834355749804793e-3, -0.341121829202314175114001e-4, 0.352995825274807520842343e-5, -0.363406184381774716643650e-6, 0.372229526012341062590015e-7, -0.379364878249523621876074e-8, 0.384735488575639572623668e-9, -0.388290050968916094626051e-10, 0.390003325917592380115783e-11, -0.389876147812208153377080e-12, 0.387934834839215670265151e-13, -0.384230036485146366813934e-14, 0.378835071251161548265558e-15, 0.290967733007319086389759e-1, -0.298556111916576121737910e-2, 0.304770397225765571278607e-3, -0.309541399850768841713903e-4, 0.312821089335649521809934e-5, -0.314583397231753700811505e-6, 0.314824345083794713386598e-7, -0.313561615787721666561815e-8, 0.310833635959190753985094e-9, -0.306698221694692246725122e-10, 0.301230841554661010546754e-11, -0.294522557960839208518663e-12, 0.286677715825897672211384e-13, -0.277811452498694880218947e-14, 0.268047105092742720807230e-15, 0.286376038148188107083594e-1, -0.289254070657916538406300e-2, 0.290707153199412470803501e-3, -0.290734726759849528451221e-4, 0.289356579447491178884808e-5, -0.286612063973937932781129e-6, 0.282558753502275012472521e-7, 
 -0.277270697024445311084857e-8, 0.270836376141295565923431e-9, -0.263356441836472481284136e-10, 0.254941302688521661220478e-11, -0.245708634307259187471724e-12, 0.235780878409829333408417e-13, -0.225282796738857839015971e-14, 0.214339139464943524401904e-15};

        private static double intExpATSqrBT(
          double alpha,
          double beta,
          double y1,
          double y2)
        {
            double bOver2a;
            double aSqr;
            double aSqry1b;
            double aSqry2b;
            double ay1bOver2a;
            double ay2bOver2a;
            double erfAB1;
            double erfAB2;
            double expAB1;
            double expAB2;
            double intgr_prev_prev;
            bOver2a = 0.5e0 * beta / alpha;
            aSqr = alpha * alpha;
            aSqry1b = y1 * (aSqr * y1 + beta);
            aSqry2b = y2 * (aSqr * y2 + beta);
            ay1bOver2a = alpha * y1 + bOver2a;
            ay2bOver2a = alpha * y2 + bOver2a;
            erfAB1 = normaldistr.erf(ay1bOver2a);
            erfAB2 = normaldistr.erf(ay2bOver2a);
            expAB1 = Math.Exp(-aSqry1b);
            expAB2 = Math.Exp(-aSqry2b);
            intgr_prev_prev = 0.5e0 * PiSqrt / alpha * Math.Exp(bOver2a * bOver2a) * (erfAB2 - erfAB1);
            return 0.5e0 / aSqr * (expAB1 - expAB2 - beta * intgr_prev_prev);
        }

        private static bool isNoNeedOfRecur(
          double b,
          double m,
          double b2OverA,
          double aSqr,
          double ba2,
          double z1,
          double z2,
          double eps0_025)
        {
            if (0.0e0 <= b)
                if (0.0e0 <= z1)
                    if (m * intExpATSqrBT(1, 0, z1, z2) <= eps0_025)
                        return true;
                    else
                        return false;
                else if (z1 < 0.0e0 && 0.0e0 <= z2 && -b2OverA <= z1)
                    if (m * (intExpATSqrBT(aSqr, -ba2, 0, -z1) + intExpATSqrBT(1, 0, 0, z2)) <= eps0_025)
                        return true;
                    else
                        return false;
                else if (z1 < 0.0e0 && 0.0e0 <= z2 && z1 < -b2OverA)
                    if (m * (intExpATSqrBT(1, 0, b2OverA, -z1) + intExpATSqrBT(aSqr, -ba2, 0, b2OverA) + intExpATSqrBT(1, 0, 0, z2)) <= eps0_025)
                        return true;
                    else
                        return false;
                else if (-b2OverA <= z1)
                    if (m * intExpATSqrBT(aSqr, -ba2, -z2, -z1) <= eps0_025)
                        return true;
                    else
                        return false;
                else if (z1 < -b2OverA && -b2OverA < z2)
                    if (m * (intExpATSqrBT(1, 0, b2OverA, -z1) + intExpATSqrBT(aSqr, -ba2, -z2, b2OverA)) <= eps0_025)
                        return true;
                    else
                        return false;
                else
                    if (m * intExpATSqrBT(1, 0, -z2, -z1) <= eps0_025)
                        return true;
                    else
                        return false;
            else
                if (z2 <= 0.0e0)
                    if (m * intExpATSqrBT(1, 0, -z2, -z1) <= eps0_025)
                        return true;
                    else
                        return false;
                else if (z1 < 0.0e0 && 0.0e0 < z2 && z2 <= b2OverA)
                    if (m * (intExpATSqrBT(1, 0, 0, -z1) + intExpATSqrBT(aSqr, ba2, 0, z2)) <= eps0_025)
                        return true;
                    else
                        return false;
                else if (z1 < 0.0e0 && b2OverA < z2)
                    if (m * (intExpATSqrBT(1, 0, 0, -z1) + intExpATSqrBT(aSqr, ba2, 0, b2OverA) + intExpATSqrBT(1, 0, b2OverA, z2)) <= eps0_025)
                        return true;
                    else
                        return false;
                else if (z2 <= b2OverA)
                    if (m * intExpATSqrBT(aSqr, ba2, z1, z2) <= eps0_025)
                        return true;
                    else
                        return false;
                else if (z1 < b2OverA && b2OverA < z2)
                    if (m * (intExpATSqrBT(aSqr, ba2, z1, b2OverA) + intExpATSqrBT(1, 0, b2OverA, z2)) <= eps0_025)
                        return true;
                    else
                        return false;
                else
                    if (m * intExpATSqrBT(1, 0, z1, z2) <= eps0_025)
                        return true;
                    else
                        return false;
        }

        private static bool isPracticalZero(
          double a,
          double b,
          double z1,
          double z2,
          double eps0_025)
        {
            if (normaldistr.erf(Math.Max(Math.Abs(a * z1 + b), Math.Abs(a * z2 + b))) * (normaldistr.erf(z2) - normaldistr.erf(z1)) <= eps0_025)
                return true;
            else
                return false;
        }

        private static bool condEps(double v, double eps)
        {
            if (eps < 0.1e1)
                return (bool)(eps < Math.Abs(v));
            else
                return true;
        }

        private static double ErfExpIntStAfter0_5(
          int iNum,
          double low,
          double up,
          double a,
          double b,
          double alpha,
          double alpha1,
          double alpha1Sqr,
          double inv2aSqr,
          double beta,
          double m,
          double b2OverA,
          double aSqr,
          double ba2,
          double eps0_025,
          double epsilon,
          int nIter)
        {
            int i;
            int indCondEpsFail;
            int indAdd;
            bool cond1;
            bool cond2;
            bool term_prev_prev_defined;
            bool term_prev_defined;
            double y1;
            double y2;
            double z1;
            double z2;
            double sh;
            double beta1;
            double bsha;
            double expMult;
            double val;
            double eps;
            double epsSqr;
            double term_prev_prev = 0.0;
            double term_prev = 0.0;
            double term;
            double bOver2a;
            double aSqry1b;
            double aSqry2b;
            double ay1bOver2a;
            double ay2bOver2a;
            double erfAB1;
            double erfAB2;
            double expAB1;
            double expAB2;
            double intgr_prev_prev;
            double intgr_prev;
            double intgr;
            double sum;
            double powy1;
            double powy2;
            z1 = alpha * low + beta;
            z2 = alpha * up + beta;
            if (isPracticalZero(a, b, z1, z2, eps0_025))
                return 0.0e0;
            if (isNoNeedOfRecur(b, m, b2OverA, aSqr, ba2, z1, z2, eps0_025))
                return normaldistr.erf(b) * (normaldistr.erf(z2) - normaldistr.erf(z1));
            sh = shifts[iNum - 1];
            y1 = low - sh;
            y2 = up - sh;
            bsha = beta + alpha * sh;
            beta1 = 0.2e1 * sh + 0.2e1 * alpha * bsha;
            expMult = Math.Exp(-sh * sh - bsha * bsha);
            eps = 0.25e-3 * epsilon / expMult;
            epsSqr = eps * eps;
            bOver2a = 0.5e0 * beta1 / alpha1;
            aSqry1b = y1 * (alpha1Sqr * y1 + beta1);
            aSqry2b = y2 * (alpha1Sqr * y2 + beta1);
            ay1bOver2a = alpha1 * y1 + bOver2a;
            ay2bOver2a = alpha1 * y2 + bOver2a;
            erfAB1 = normaldistr.erf(ay1bOver2a);
            erfAB2 = normaldistr.erf(ay2bOver2a);
            expAB1 = Math.Exp(-aSqry1b);
            expAB2 = Math.Exp(-aSqry2b);
            cond1 = false;
            cond2 = false;
            term_prev_prev_defined = false;
            term_prev_defined = false;
            indAdd = nIter * (iNum - 1);
            intgr_prev_prev = 0.5e0 * PiSqrt / alpha1 * Math.Exp(bOver2a * bOver2a) * (erfAB2 - erfAB1);
            val = coeffsErfc[indAdd] * intgr_prev_prev;
            sum = val;
            if (condEps(val, epsSqr))
            {
                term_prev_prev = val;
                term_prev_prev_defined = true;
                cond1 = (bool)(Math.Abs(term_prev_prev) <= eps);
            }
            intgr_prev = 0.5e0 / alpha1Sqr * (expAB1 - expAB2 - beta1 * intgr_prev_prev);
            val = coeffsErfc[indAdd + 1] * intgr_prev;
            if (condEps(val, epsSqr))
            {
                sum += val;
                if (term_prev_prev_defined)
                {
                    term_prev = val;
                    term_prev_defined = true;
                }
                else
                {
                    term_prev_prev = val;
                    term_prev_prev_defined = true;
                    cond1 = (bool)(Math.Abs(term_prev_prev) <= eps);
                }
            }
            powy1 = 0.1e1;
            powy2 = 0.1e1;
            indCondEpsFail = 0;
            for (i = 2; i <= nIter - 1; i++)
            {
                if (3 <= indCondEpsFail)
                    return normaldistr.erf(z2) - normaldistr.erf(z1) - 0.4e1 * alpha / PiSqrt * expMult * sum;
                powy1 *= y1;
                powy2 *= y2;
                intgr = inv2aSqr * ((double)(i - 1) * intgr_prev_prev - beta1 * intgr_prev + expAB1 * powy1 - expAB2 * powy2);
                intgr_prev_prev = intgr_prev;
                intgr_prev = intgr;
                val = coeffsErfc[indAdd + i] * intgr;
                if (condEps(val, epsSqr))
                {
                    indCondEpsFail = 0;
                    sum += val;
                    if (term_prev_prev_defined)
                        if (term_prev_defined)
                        {
                            term = val;
                            cond2 = (bool)(Math.Abs(term) / 0.3e1 + Math.Abs(term_prev) / 0.3e1 + Math.Abs(term_prev_prev) / 0.3e1 <= eps);
                            if (cond1 && cond2)
                                return normaldistr.erf(z2) - normaldistr.erf(z1) - 0.4e1 * alpha / PiSqrt * expMult * sum;
                            term_prev_prev = term_prev;
                            term_prev = term;
                            cond1 = (bool)(Math.Abs(term_prev_prev) <= eps);
                        }
                        else
                        {
                            term_prev = val;
                            term_prev_defined = true;
                        }
                    else
                    {
                        term_prev_prev = val;
                        term_prev_prev_defined = true;
                        cond1 = (bool)(Math.Abs(term_prev_prev) <= eps);
                    }
                }
                else
                    indCondEpsFail++;
            }
            return normaldistr.erf(z2) - normaldistr.erf(z1) - 0.4e1 * alpha / PiSqrt * expMult * sum;
        }

        private static double ErfExpIntStBefore0_5(
          double low,
          double up,
          double a,
          double b,
          double alpha,
          double beta,
          double m,
          double b2OverA,
          double aSqr,
          double ba2,
          double eps0_025)
        {
            int i;
            int i2;
            double z1;
            double z2;
            double alphaSqr;
            double beta1;
            double inv2aSqr;
            double aSqry1b;
            double aSqry2b;
            double ay1bOver2a;
            double ay2bOver2a;
            double erfAB1;
            double erfAB2;
            double expAB1;
            double expAB2;
            double intgr_even;
            double intgr_odd;
            double sum;
            double pLow;
            double pUp;
            z1 = alpha * low + beta;
            z2 = alpha * up + beta;
            if (isPracticalZero(a, b, z1, z2, eps0_025))
                return 0.0e0;
            if (isNoNeedOfRecur(b, m, b2OverA, aSqr, ba2, z1, z2, eps0_025))
                return normaldistr.erf(b) * (normaldistr.erf(z2) - normaldistr.erf(z1));
            beta1 = 0.2e1 * alpha * beta;
            alphaSqr = alpha * alpha;
            inv2aSqr = 0.5e0 / alphaSqr;
            aSqry1b = low * (alphaSqr * low + beta1);
            aSqry2b = up * (alphaSqr * up + beta1);
            ay1bOver2a = alpha * low + beta;
            ay2bOver2a = alpha * up + beta;
            erfAB1 = normaldistr.erf(ay1bOver2a);
            erfAB2 = normaldistr.erf(ay2bOver2a);
            expAB1 = Math.Exp(-aSqry1b);
            expAB2 = Math.Exp(-aSqry2b);
            intgr_even = 0.5e0 * PiSqrt * a * Math.Exp(beta * beta) * (erfAB2 - erfAB1);
            intgr_odd = inv2aSqr * (expAB1 - expAB2 - beta1 * intgr_even);
            sum = intgr_odd;
            pLow = 0.1e1;
            pUp = 0.1e1;
            for (i = 1; i <= 11; i++)
            {
                i2 = 2 * i;
                pLow *= low;
                pUp *= up;
                intgr_even = inv2aSqr * ((double)(i2 - 1) * intgr_even - beta1 * intgr_odd + expAB1 * pLow - expAB2 * pUp);
                pLow *= low;
                pUp *= up;
                intgr_odd = inv2aSqr * ((double)i2 * intgr_odd - beta1 * intgr_even + expAB1 * pLow - expAB2 * pUp);
                sum += coeffsErf[i - 1] * intgr_odd;
            }
            return 0.4e1 / Math.PI * alpha * sum * Math.Exp(-beta * beta);
        }

        private static int getPos(double y)
        {
            if (y <= 0.5e0)
                return 0;
            if (0.5e0 < y && y <= 0.7375e0)
                return 1;
            if (0.7375e0 < y && y <= 0.975e0)
                return 2;
            if (0.975e0 < y && y <= 0.12125e1)
                return 3;
            if (0.12125e1 < y && y <= 0.145e1)
                return 4;
            if (0.145e1 < y && y <= 0.16875e1)
                return 5;
            if (0.16875e1 < y && y <= 0.1925e1)
                return 6;
            if (0.1925e1 < y && y <= 0.21625e1)
                return 7;
            if (0.21625e1 < y && y <= 0.24e1)
                return 8;
            if (0.24e1 < y && y <= 0.26375e1)
                return 9;
            if (0.26375e1 < y && y <= 0.2875e1)
                return 10;
            if (0.2875e1 < y && y <= 0.31125e1)
                return 11;
            if (0.31125e1 < y && y <= 0.335e1)
                return 12;
            if (0.335e1 < y && y <= 0.35875e1)
                return 13;
            if (0.35875e1 < y && y <= 0.3825e1)
                return 14;
            if (0.3825e1 < y && y <= 0.40625e1)
                return 15;
            if (0.40625e1 < y && y <= 0.43e1)
                return 16;
            if (0.43e1 < y && y <= 0.45375e1)
                return 17;
            if (0.45375e1 < y && y <= 0.4775e1)
                return 18;
            if (0.4775e1 < y && y <= 0.50125e1)
                return 19;
            if (0.50125e1 < y && y <= 0.525e1)
                return 20;
            if (0.525e1 < y && y <= 0.54875e1)
                return 21;
            if (0.54875e1 < y && y <= 0.5725e1)
                return 22;
            if (0.5725e1 < y && y <= 0.59625e1)
                return 23;
            if (0.59625e1 < y && y <= 0.62e1)
                return 24;
            if (0.62e1 < y && y <= 0.64375e1)
                return 25;
            if (0.64375e1 < y && y <= 0.6675e1)
                return 26;
            if (0.6675e1 < y && y <= 0.69125e1)
                return 27;
            if (0.69125e1 < y && y <= 0.715e1)
                return 28;
            if (0.715e1 < y && y <= 0.73875e1)
                return 29;
            if (0.73875e1 < y && y <= 0.7625e1)
                return 30;
            if (0.7625e1 < y && y <= 0.78625e1)
                return 31;
            if (0.78625e1 < y && y <= 0.81e1)
                return 32;
            if (0.81e1 < y && y <= 0.83375e1)
                return 33;
            if (0.83375e1 < y && y <= 0.8575e1)
                return 34;
            if (0.8575e1 < y && y <= 0.88125e1)
                return 35;
            if (0.88125e1 < y && y <= 0.905e1)
                return 36;
            if (0.905e1 < y && y <= 0.92875e1)
                return 37;
            if (0.92875e1 < y && y <= 0.9525e1)
                return 38;
            if (0.9525e1 < y && y <= 0.97625e1)
                return 39;
            if (0.97625e1 < y && y <= 0.100e2)
                return 40;
            return 41;
        }

        private static double ErfExpIntPosBounds(
          double sourceLow,
          double sourceUp,
          double low,
          double up_,
          double a,
          double b,
          double epsilon)
        {
            double m;
            double b2OverA;
            double aSqr;
            double ba2;
            double alpha;
            double beta;
            double eps0_025;
            double up;
            double alpha1;
            double alpha1Sqr;
            double inv2aSqr;
            double sum;
            int lowPos;
            int upPos;
            int i;
            m = 0.4e1 / Math.PI * Math.Exp(-b * b) * a;
            b2OverA = 0.2e1 * Math.Abs(b) / a;
            aSqr = Math.Sqrt(0.1e1 + a * a);
            ba2 = 0.2e1 * a * b;
            alpha = 0.1e1 / a;
            beta = -b * alpha;
            eps0_025 = 0.25e-1 * epsilon;
            if (isPracticalZero(a, b, sourceLow, sourceUp, epsilon))
                return 0.0e0;
            if (isNoNeedOfRecur(b, m, b2OverA, aSqr, ba2, sourceLow, sourceUp, epsilon))
                return normaldistr.erf(b) * (normaldistr.erf(sourceUp) - normaldistr.erf(sourceLow));
            lowPos = getPos(low);
            upPos = getPos(up_);
            if (nBounds < upPos)
            {
                upPos = nBounds;
                up = 0.100e2;
                if (nBounds < lowPos)
                    return normaldistr.erf(sourceUp) - normaldistr.erf(sourceLow);
                else
                    sum = normaldistr.erf(sourceUp) - normaldistr.erf(alpha * up + beta);
            }
            else
            {
                up = up_;
                sum = 0.0e0;
            }
            if (upPos == 0)
                return ErfExpIntStBefore0_5(low, up, a, b, alpha, beta, m, b2OverA, aSqr, ba2, eps0_025);
            else
            {
                alpha1 = Math.Sqrt(0.1e1 + alpha * alpha);
                alpha1Sqr = alpha1 * alpha1;
                inv2aSqr = 0.5e0 / alpha1Sqr;
                if (lowPos == 0)
                    sum += ErfExpIntStBefore0_5(low, 0.5e0, a, b, alpha, beta, m, b2OverA, aSqr, ba2, eps0_025);
                else
                    if (lowPos == upPos)
                    {
                        sum += ErfExpIntStAfter0_5(lowPos, low, up, a, b, alpha, alpha1, alpha1Sqr, inv2aSqr, beta, m, b2OverA, aSqr, ba2, eps0_025, epsilon, nIter);
                        return sum;
                    }
                    else
                        sum += ErfExpIntStAfter0_5(lowPos, low, bounds[lowPos], a, b, alpha, alpha1, alpha1Sqr, inv2aSqr, beta, m, b2OverA, aSqr, ba2, eps0_025, epsilon, nIter);
                for (i = lowPos + 1; i <= upPos - 1; i++)
                    sum += ErfExpIntStAfter0_5(i, bounds[i - 1], bounds[i], a, b, alpha, alpha1, alpha1Sqr, inv2aSqr, beta, m, b2OverA, aSqr, ba2, eps0_025, epsilon, nIter);
                sum += ErfExpIntStAfter0_5(upPos, bounds[upPos - 1], up, a, b, alpha, alpha1, alpha1Sqr, inv2aSqr, beta, m, b2OverA, aSqr, ba2, eps0_025, epsilon, nIter);
                return sum;
            }
        }

        private static double ErfExpIntSpecial(double a, double b, double x)
        {
            double low;
            double up;
            double q;
            low = -a * xT + b;
            up = a * x + b;
            if (up <= 0.0e0)
                return -ErfExpIntPosBounds(-x, xT, -up, -low, a, -b, epsHalf);
            else if (low < 0.0e0)
            {
                q = b / a;
                return -ErfExpIntPosBounds(q, xT, 0, -low, a, -b, epsQuarter) + ErfExpIntPosBounds(-q, x, 0, up, a, b, epsQuarter);
            }
            else
                return ErfExpIntPosBounds(-xT, x, low, up, a, b, epsHalf);
        }
    }
}
